Upper limit on the ultra-high-energy photon flux from AGASA and Yakutsk data 
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We present the interpretation of the muon and scintillation signals of ultra-high-energy air showers 
observed by AGASA and Yakutsk extensive air shower array experiments. We consider case-by-case 
ten highest energy events with known muon content and conclude that at the 95% confidence level 
(C.L.) none of them was induced by a primary photon. Taking into account statistical fluctuations 
and differences in the energy estimation of proton and photon primaries, we derive an upper limit 
of 36% at 95% C.L. on the fraction of primary photons in the cosmic-ray flux above 10 20 eV. This 
result disfavors the Z-burst and superheavy dark-matter solutions to the GZK-cutoff problem. 



PACS numbers: 98.70.Sa, 96.40.De, 96.40.Pq 

I. INTRODUCTION 

One of the most intriguing puzzles in astroparticle 
physics is the observation of air showers initiated by 
particles with energies beyond the cutoff predicted by 
Greisen and by Zatsepin and Kuzmin Compared 
to lower energies, the energy losses of protons increase 
sharply at « 5 x 10 19 eV since pion production on cosmic 
microwave background photons reduces the proton mean 
free path by more than two orders of magnitude. This 
effect is even stronger for heavier nuclei, while photons 
are absorbed due to pair production on the radio back- 
ground with the mean free path of a few Mpc. Thus, 
the cosmic-ray (CR) energy spectrum should dramati- 
cally steepen at « 7 x 10 19 eV for any homogeneous dis- 
tribution of CR sources. Despite the contradictions in 
the shape of the spectrum, the existence of air showers 
with energies in excess of 10 20 eV is firmly established 
by several independent experiments using different tech- 
niques (Volcano Ranch 0, Fly's Eye Q, Yakutsk 0, 
AGASA 0, HiRes and Pierre Auger 13 experiments). 
Some explanations for these showers, like the Z-burst or 
top-down models, predict a significant fraction of pho- 
tons above typically 8 x 10 19 eV (for reviews see, e.g., 
Refs. Q). Indications for the presence of neutral par- 
ticles at lower energies were found in Refs. 0- Thus, 
the determination of the photon fraction in the CR flux 
is of crucial importance, and the aim of this work is to 
derive a stringent limit on this fraction in the integral 
CR flux above 10 20 eV. To this end, we compare the re- 
ported information on signals measured by scintillation 
and by muon detectors for observed showers with those 
expected by air shower simulations. We focus on the sur- 
face detector signal density at 600 meters S^OO) (known 
as charged particle density) and the muon density at 1000 
m, /^(lOOO), which are used in experiments as primary 
energy and primary mass estimators, respectively. 

We study individual events of AGASA 0] and of the 
Yakutsk extensive air shower array (Yakutsk in what fol- 



lows) Q with reconstructed energies above 8 x 10 19 eV 
and measured muon content. We reject the hypothe- 
sis that any of showers considered was initiated by a 
photon primary at the 95% confidence level (C.L.). We 
then derive as our main result an upper limit of 36% 
(at 95% C.L.) on the fraction e 7 of primary photons 
with original energies above 10 20 eV (the difference be- 
tween original and reconstructed energies is discussed in 
Sec.|nJ. 

The rest of the paper is organized as follows. In SecHTI 
we discuss the experimental data set which we use for 
our study. In Sec. IIIII the details of the simulation of 
the artificial shower libraries and comparison of the sim- 
ulated and real data are given. This section contains the 
description of our method and the main results. We dis- 
cuss how robust these results are with respect to changes 
in assumptions, to analysis procedure, and to variations 
in the experimental data, in Sec. lIVI In Sec.|Vl we discuss 
the differences between our approach and previous stud- 
ies, which allowed us to put a significantly more strin- 
gent limit on the gamma-ray fraction. Our conclusions 
are briefly summarized in Sec. IVII 



II. EXPERIMENTAL DATA 

AGASA was operating from 1990 to 2003 and consisted 
of 111 surface scintillation detectors (covering an area of 
about 100 km 2 ) and 27 muon detectors. The areas of the 
AGASA muon detectors varied between 2.8 and 20 m 2 . 
The detectors consisted of 14-20 proportional counters 
aligned under a shield of either 30 cm of iron or 1 m of 
concrete and were placed below or close to scintillation 
detectors. The threshold energy was 0.5 GeV/ cos# M for 
muons with zenith angle 9^ 11]. During 14 years of op- 
eration, AGASA had observed 11 events with reported 
energies above 10 20 eV and zenith angles 9 < 45° !5l[T2l. 
Among them, six events had p M (1000) determined \T\\ . 

Yakutsk is observing CRs of highest energies since 
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1973, with detectors in various configurations. With 
6 < 60°, it has observed three events above 10 20 eV, 
all with measured muon content. Before 1978, only one 
muon detector with the area of 8 m 2 and threshold en- 
ergy 0.7 GeV/ cos6* M was in operation. Later, it has been 
replaced by six detectors with areas up to 36 m 2 and the 
threshold energy of 1.0 GeV/ cos 9^ |l3| . 

In our study, we combine the AGASA and Yakutsk 
datasets, motivated by the following. First, both datasets 
are obtained from surface array experiments operated 
with similar plastic scintillation detectors. Second, the 
energy estimation procedures of the two experiments 
are compatible, within the reported systematic errors at 
~ 10 20 eV, if differences in the observational conditions 
are taken into account 0]. Finally, the values of the 
CR flux at 10 20 eV reported by the two experiments are 
consistent within their la errors. 

The shower energy estimated by an experiment (here- 
after denoted as E est ) is in general different from the 
true primary energy (denoted as Eo) because of natu- 
ral shower fluctuations, etc. Moreover, the energy es- 
timation algorithms used by surface-array experiments 
normally assume that the primary is a proton. While 
the estimated energy for nuclei depends only weakly on 
their mass number, the difference between photons and 
hadrons is significant. For photons, the effects of geo- 
magnetic field [L5j result in directional dependence of the 
energy reconstruction. Thus, the event energy reported 
by the experiment should be treated with care when we 
allow the primary to be a photon. In this study we in- 
clude events with E cst > 8 x 10 19 eV because of pos- 
sible energy underestimation for photon-induced show- 
ers; these events contribute to the final limit, derived for 
E Q > 10 20 eV, with different weights. 

For AGASA, we use the events given in Ref. [B[ that 
pass the "cut B" defined in Ref. HU, that is having 
at least one 46] muon detector hit between 800 m and 
1600 m from the shower axis. The /^(lOOO) of the in- 
dividual events can be read off from Fig. 2 of Ref. [IT| . 
Yakutsk muon detectors have larger area and are more 
sensitive both to weak signals far from the core and to 
strong signals for which AGASA detectors might become 
saturated. This allowed the Yakutsk collaboration to re- 
lax the cuts, as compared to AGASA, and to obtain re- 
liable values of ,0^(1000) using detectors between 400 m 
and 2000 m from the shower axis 0, 0|. Providing 
these cuts, six AGASA and four Yakutsk events entered 
the dataset in our study (see Table [I] for the event de- 
tails). 



III. SIMULATIONS AND RESULTS 

In order to interpret the data, for each of the ten 
events, we generated a shower library containing 1000 
showers induced by primary photons 47j . Thrown ener- 
gies Eq of the simulated showers were randomly selected 
(see below the discussion of the initial spectra) between 



5 x 10 19 eV and 5 x 10 20 eV to take into account possible 
deviations of E est from Eq . The arrival directions of the 
simulated showers were the same as those of the corre- 
sponding real events. The simulations were performed 
with CORSIKA v6.204 [3, choosing QGS JET 01c 
as high-energy and FLUKA 2003.1b [2(| as low-energy 
hadronic interaction model. Electromagnetic shower- 
ing was implemented with EGS4 |2l| incorporated into 
CORSIKA. Possible interactions of the primary pho- 
tons with the geomagnetic field were simulated with the 
PRES HOWE R option of CORSIKA 22]. As discussed 
in Sec. lIVBl this choice of the interaction models results 
in a conservative limit on gamma-ray primaries. As sug- 
gested in Ref. .23]], all simulations were performed with 
thinning level 10 -5 , maximal weight 10 6 for electrons and 
photons, and 10 for hadrons. 

For each simulated shower, we determined ,5(600) and 
p M (1000). For the calculation of 5(600), we used the de- 
tector response functions from Refs. 0, |2j| . For a given 
arrival direction, there is one-to-one correspondence be- 
tween 5(600) and the quantity called estimated energy, 
-B cs t • The relation is determined b y th e standard analysis 
procedure of the two experiments [Hj, |22| • This allows us 
to select simulated showers compatible with the observed 
ones by the signal density. The quantity 5(600) is recon- 
structed not precisely. In terms of estimated energy, for 
AGASA events, the reconstructed energies are are dis- 
tributed with a Gaussian in log (E^t/ E rpr )\ the standard 
deviation of E est is a s=s 25% [lj]. For Yakutsk events, 
the corresponding a has been determined event-by-event 
and is typically 30-45% j2?|. To each simulated shower, 
we assigned a weight wi proportional to this Gaussian 
probability distribution in log E cst centered at the ob- 
served energy i? roc = E ^ s . Additionally, each simulated 
shower was weighted with u>2 to reproduce the thrown 
energy spectrum oc Eq 2 which is typically predicted by 
non-acceleration scenarios (see Sec. IIV 01 for a discussion 
of the variations of the spectral index) . For each of the 
ten observed events, we obtained a distribution of muon 
densities p M (1000) representing photon-induced showers 
compatible with the observed ones by 5(600) and arrival 
directions. To this end, we calculated p M (1000) for each 
simulated shower by making use of the same muon lat- 
eral distribution function as used in the analysis of real 
data 0, 0] . To take into account possible experimen- 
tal errors in the determination of the muon density, we 
replaced each simulated /o M (1000) by a distribution repre- 
senting possible statistical errors (50% and 25% Gaussian 
for AGASA cut B 29] and Yakutsk [l7|, respectively). 
The distribution of the simulated muon densities is the 
sum of these Gaussians weighted by w\W2- 

A typical distribution of simulated p M (1000) is given in 
Fig-ffl f° r gamma- and proton- induced simulated showers 
compatible with the event 3 by 5(600) and the arrival di- 
rection. We will see below that this particular event has 
the largest probability of gamma interpretation among 
all ten events in the data set; still the proton interpreta- 
tion looks perfect for it. This is the case for all events 
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TABLE I: Description of the individual events used in this work. Columns: (1), event number; (2), experiment; (3), date of the 
event detection (in the format dd.mm.yyyy); (4), the reported energy assuming a hadronic primary (in units of 10 20 eV); (5), 
the zenith angle (in degrees); (6) the azimuth angle (in degrees, (f> — corresponds to a particle coming from the South, <j) = 90° 
- from the West); (7) number of muon detectors used to reconstruct muon density; (8) muon density at 1000 m from the shower 
axis (in units of m -2 ); (9), probability that this event was initiated by a photon with E > 10 20 eV; (10), probability that this 
event was initiated by a non-photon with E > 10 20 eV, assuming correct energy determination. The sum p^ + p^ gives the 
weight of this event in the final limit on e 7 . The probability that the primary had the energy E < 10 20 eV is 1 — p^' — p^ ■ 
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FIG. 1: Weighted distributions of muon density p M (1000) for 
the simulated events compatible with the event 3 by 5(600) 
and the arrival direction. Units in the vertical axis are ar- 
bitrary, p M (1000) is measured in m -2 . The thin dark line 
corresponds to primary photons; it is the distribution used 
for our analysis. The thick grey line is the distribution ob- 
tained in the same way but for 500 proton-induced showers. 
The arrow indicates the observed value of p M (1000) for the 
event 3. The disributions include 50% Gaussian error of the 
detector. 



except event 7, which has too high p M (1000) for a proton; 
possible nature of its primary particle will be discussed 
elsewhere. 

To estimate the allowed fraction e 7 of primary photons 
among CRs with Eo > 10 20 eV, we compare, for each ob- 
served event, two possibilities: (i) that it was initiated by 
a photon primary with Eq > 10 20 eV and (ii) that it was 
initiated by any other primary with Eq > 10 20 eV for 
which the experimental energy estimation works prop- 
erly. 



Let us consider the ith observed event. Denote by 
M the weighted number of showers contributed to the 
p M (1000) distribution for the simulated photon-induced 
showers compatible with the ith event by arrival direction 
and 5(600) (throughout this paragraph, the weighted 
number is the sum of corresponding weights, that is 
M is the sum of weights of all 1000 showers simu- 
lated for the ith event). Some of the simulated show- 
ers contributed to the part of the distribution for which 
/0 M (1000) > ^(1000), where p^(1000) is the observed 
value for this event. The weighted number of these show- 
ers is M'. Some part I of this M' corresponds to showers 
with E > 10 20 eV, the rest (M' - I) to E Q < 10 20 eV. 

The probability p^ of case (i) is pf* = l/M, while the 
probability that the event is consistent with a photon 
of E < 10 20 eV is pf ] = (M' - l)/M. Moreover, 
the probability that the event is described by any other 
primary is 1 — p^ — p'^ 1 = 1 — M'/M. We assume 
that the experimental energy estimation works well for 
non-photon primaries and determine the fraction £ of 
events with E > 10 20 eV simply from the Gaussian 
log(i?ost) distribution, so the probability of the case (ii) 
is p^ = £(1 — M'/M). The values of p± 2 are presented 

in Table [I] Note that py + p% < 1 because of a non- 
zero probability that a simulated shower is initiated by 
a primary with Eq < 10 20 eV. This happens especially 
for events with reported energies close to 10 20 eV and 
reduces considerably the effective number of events con- 
tributing to the limit on e 7 : since we are interested in 
the limit for Eo > 10 20 eV only, each event contributes 
to the result with the weight (p^ +i4 )■ Inspection of 
Table [I] demonstrates that the total effective number of 
events with Eo > 10 20 eV (the sum of pf 1 and p% over 
all ten events) is 7.67. 
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If the ith primary particle was a photon with Eq > 
10 20 eV with the probability p± and a non-photon with 
Eq > 10 20 eV with the probability p% , one can eas- 
ily calculate the probability V{ni,n 2 ) to have ri\ pho- 
tons and n 2 non-photons in the set of N = 10 observed 
events (0 < n\ + n 2 < N, the rest N — rii — n 2 events 
have Eq < 10 20 eV). From the set of N events, one 
should take all possible non-overlapping subsets of n\ 
and n 2 events and sum up probabilities of these real- 
isations (since p\\ ^ V\\i these probabilities are dif- 
ferent for different realisations with the same n\ and 
ri2). Now, suppose that the fraction of the primary pho- 
tons at E > 10 20 eV is e 7 . Then, the probability to 
have n\ photons and n 2 non-photons at Eq > 10 20 eV 
is e™ 1 (1 — e 7 )" 2 , and the probability that the observed 
muon densities were obtained with a given e 7 is 

N 

n\ ,^2—0 

(cf. Ref. [3(j for a particular case ri\ + n 2 = N; note that 
the combinatorial factor is included in the definition of 
V(n%, n 2 )). The cases n\ + n 2 < N reflect the possibility 
that some of the N events correspond to primaries with 
Eq < 10 20 eV. In our case, the probability V(e 7 ) is a 
monotonically decreasing function of e 7 . Thus the up- 
per limit on e 7 at the confidence level a' is obtained by 
solving the equation 7 ? (e 7 ) = 1 — a'. For our dataset, the 
95% C.L. upper limit on the photon fraction is e 7 < 0.33. 
The limit on e 7 is rather weak compared to the individ- 
H) 

ual values of p\ because of the small number of observed 
events. 

However, some of the photon-induced showers may es- 
cape from our study because they may not pass the muon 
measurement quality cuts or their estimated energy is be- 
low 8 x 10 19 eV. Possible reasons for an underestimation 
of the energy may be either the LPM effect [3l| or sub- 
stantial attenuation of gamma-induced showers at large 
zenith angles. To estimate the fraction of these "lost" 
events, we have simulated 1000 gamma-induced showers 
for each experiment with arrival directions distributed 
according to the experimental acceptance. We find that 
the fraction of the "lost" events is ~ 3.5% for AGASA 
and ~ 15% for Yakutsk. The account of these fractions, 
weighted with the relative exposures of both experiments, 
results in the final upper limit, 

e 7 < 36% (95% C.L.). 

In Fig. [3 we present our limit on e 7 (AY) together 
with previously published limits on the same quantity. 
Also, typical theoretical predictions are shown for the 
superheavy dark-matter, topological-defect and Z-burst 
models. Our limit on e 7 is currently the strongest one 
at Eq > 10 20 eV. It disfavours some of the theoretical 
models such as the Z-burst and superheavy dark-matter 
scenarios. 
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FIG. 2: Limits (95% C.L.) on the fraction e 7 of photons in 
the integral CR flux versus energy. The result of the present 
work (AY) is shown together with limits previously given in 
Refs. |I3 (HP), E3 (A), (RH) and H3 (PA). Also shown 
are predictions for the superheavy dark-matter (thick line) 
and topological-defect (necklaces, between dotted lines) mod- 
els |34j) and for the Z-burst model (shaded area) |3Sll . 



IV. ROBUSTNESS OF THE RESULTS 

In this section, we discuss systematic uncertainties of 
our limit that are related to the air shower simulations, 
to the data interpretation and to selection cuts. 



A. Systematic uncertainty in the ,5(600) and energy 
determination 

The systematic uncertainty in the absolute energy 
determination is 18% and 30% for AGASA [l(]] and 
Yakutsk respectively. These systematic errors origi- 
nate from two quite different sources: (a) the measure- 
ment of 5(600) and (b) the relation between S*(600) and 

(i) 

primary energy. The probabilities p\ that a particu- 
lar event may allow for a gamma-ray interpretation are 
not at all sensitive to the S'(600)-to-energy conversion 
because we select simulated events by 5(600) and not 
by energy. These probabilities may be affected by rela- 
tive systematics in determination of p (tt (1000) and 5(600). 

On the other hand, in the calculation of p% we assumed 
that the experimental energy determination is correct for 
non-photon primaries; the values of p 2 ^ and the effec- 
tive number of events contributing to the limit on e 7 at 
Eq > 10 20 eV would change if the energies are systemati- 
cally shifted. In our case (all pf 1 ~ 0), the reported value 
of e 7 would be applicable to the shifted energy range in 
that case. 

Thus, the 95% C.L. conclusion that none of the ten 
events considered here was initiated by a photon is ro- 
bust with respect to any changes in the 5(600)-to-energy 
conversion. As for the limit on e 7 we report, instead of 
Eq > 10 20 eV, it would be applicable to a different en- 
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ergy range if all experimental energies are systematically 
shifted. One should note that theoretical predictions, e.g. 
the curves shown in Fig. would also change because 
they are normalised to the observed AGASA spectrum. 



B. Interaction models and simulation codes 

Our simulations were performed entirely in the COR- 
SIKA framework, and any change in the interaction mod- 
els or simulation codes, which affects either 5(600) or 
p M (1000), may affect our limit. We have studied the 
model dependence of our results by comparing differ- 
ent low- and high-energy hadronic interaction models 
(GHE1SHA HI versus FLUKA, SIBYLL 2.1 ^ ver- 
sus QGSJET). Our result is quite stable with respect to 
these changes. In all cases, individual values of are 
always close to zero, thus the limit on e 7 is not affected. 
The change of the low energy model does not at all af- 
fect the reported values. In use of SIBYLL compared 
with QGSJET, ^^(1000) is - 20% smaller for photon- 
induced showers. While 5(600) is almost unchanged, 
events in our dataset are better explained by showers 
initiated by heavier nuclei and the probability of photon- 
induced showers is even smaller. A similar effect is ex- 
pected for the coming interaction model QGSJET II [38|. 

We also performed simulations with the help of the 
hybrid code which reproduced the CORSIKA re- 
sults with high accuracy. Another popular simulation 
code, AIRES lH, differs from CORSIKA mainly in the 
low-energy hadronic interaction model (which is fixed 
in AIRES to be the Hillas splitting algorithm), hence 
we hope that simulations with AIRES would not signif- 
icantly affect our results. Comparison with AIRES will 
be presented elsewhere. 

The values presented here were obtained for the stan- 
dard parameterization of the photo-nuclear cross section 
given by the Particle Data Group (implemented as 
default in CORSIKA). The muon content of gamma- 
induced showers is in principle sensitive to the extrap- 
olation of the photonuclear cross section to high ener- 
gies. The hybrid code jl^ allows for easy variations of 
the cross section; we checked that the results are stable 
for various reasonable extrapolations, in agreement with 
Ref. [12. 



C. Primary energy spectrum 

For our limit, we used the primary photon spectrum 
E^ a for a = 2. While the individual probabilities p±\ are 
not affected by the change of the spectral index a because 
the simulated events are selected by 5(600) anyway, the 
value of a changes the fraction of "lost" photons and, 
correspondingly, the final limit on e 7 . Variations of 1 < 
a < 3 result in the photon fraction limits between 36% 
and 37% (95% C.L.). 



D. Width of the distribution 

Clearly, the rare probabilities of high values of 
p M (1000) in the tail of the distribution for primary pho- 
tons depend on the width of this distribution. The fol- 
lowing sources contribute to this width: 

• variations of the primary energy compatible with 
the observed 5(600) (larger energy correspond to 
larger muon number and ,0^(1000)); 

• physical shower-to-shower fluctuations in muon 
density for a given energy (dominated by fluc- 
tuations in the first few interactions, including 
preshowering in the geomagnetic field); 

• artificial fluctuations in 5(600) and ^(1000) due 
to thinning; 

• experimental errors in p^(lOOO) determination. 

While the first two sources are physical and are fully 
controlled by the simulation code, the variations of the 
last two may affect the results. 



1. Artificial fluctuations due to thinning 

It has been noted in Ref. |43| that the fluctuations in 
p M (1000) due to thinning may affect strongly the preci- 
sion of the composition studies. For the thinning param- 
eters we use, the relative size of these fluctuations is Q 
- 10% for ^(1000) and - 5% for 5(600). Thus with 
more precise simulations, the distributions of muon den- 
sities should become more narrow, which would reduce 
the probability of the gamma-ray interpretation of each 
of the studied events even further. 



2. Experimental errors in p M (1000) determination 

The distributions of (1000) we use accounted for the 
error in experimental determination of this quantity. The 
size of the errors was taken from the original experimental 
publications p"?l |29|. In principle, this error depends on 
the event quality and on the muon number itself, which 
is lower for simulated gamma-induced showers than for 
the observed ones. However, e.g. Ref. states that 
for the AGASA cut A (two or more muon detectors), 
the error is 40%, lower than 50% we use j2^|. Note that 
Ref. 01 discusses muon densities as low as 0.04 m~ 2 and 
even m~ 2 , much lower than ~ 1 m~ 2 typical for our 
simulated gamma-induced events. Still, we tested the 
stability of our limit by taking artificially high values of 
experimental errors in muon density: 100% for AGASA 
and 50% for Yakutsk. The limit on e 7 changes to 37% 
(95% C.L.) in that case. 
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E. Data selection cuts 

Since all events in the data set are unlikely to be ini- 
tiated by primary photons (all w 0), the limit on 
e 7 is determined by statistics only and is affected if the 
number of events is changed. Here, we discuss possible 
variations of the data set corresponding to more stringent 
quality cuts which reduce the event number and weaken 
the limit. 



1. Zenith angle 

All Yakutsk events in the data set have zenith angles 
45° < 9 < 60°, so the cut 9 < 45° imposed by AGASA 
reduces the sample to six AGASA events which results in 
the limit e 7 < 50% (95% C.L.). One should note however 
that AGASA muon detectors are not sensitive to inclined 
showers, which is not the case for Yakutsk. 



2. Core inside array 

Another cut imposed on the AGASA published dataset 
is the location of the core inside array. The event number 
7 does not satisfy this criterion; its exclusion from the 
data set results in e 7 < 40% (95% C.L.). 



3. More than one muon detector 

Reconstruction of the muon density at 1000 m from 
a single muon detector reading requires extrapolation of 
the lateral distribution function with an averaged slope. 
Though it is well-studied, the data points corresponding 
to events with a single muon detector hit might be consid- 
ered less reliable than those with two or more hits. With 
the account of the events with two or more hits only, 
we are left with seven events (four AGASA and three 
Yakutsk) which weakens the 95% C.L. limit to e 7 < 48%. 



V. COMPARISON WITH OTHER STUDIES 

Some of the previous studies used the AGASA 0, 
301 and Yakutsk muon data to limit the gamma-ray 
primaries at high energies. Our results differ from the 
previous ones not only because we join the data sets of 
the two experiments. Two major distinctive features of 
our approach allowed us to put the stringent limit: 

• both p M (1000) and 5(600) were tracked for simu- 
lated showers within framework of a single simu- 
lation code (CORSIKA in our case); 

• each event was studied individually, without aver- 
aging over arrival directions. 
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FIG. 3: Direction dependence of the reconstructed energy 
for gamma-ray primaries. Plotted is the reconstructed energy 
(determined by the AGASA method from S^OO)) versus the 
primary energy. Dark boxes: arrival direction of the event 1; 
crosses: arrival direction of the event 3; grey circles: arrival 
directions randomly distributed according to the AGASA ac- 
ceptance (0 < < 45°). Straight line represents _B roc = Eq. 
Both Eq and E ICC are measured in eV. 



In Refs. 0, 0, no conclusion was derived about e 7 
at E > 10 20 eV, and the data points corresponding to 
highest-energy events were found to be quite close to the 
gamma-ray domain. To our opinion, the main source of 
this effect is averaging over arrival directions which intro- 
duced additional fluctuations for gamma-ray primaries 
due to direction-depende nt p reshowering (see Fig. [3] for 
an illustration). In Ref. |30| which discussed the same 
six AGASA events, all simulated showers for an event 
with the observed energy -B b s had energies 1.2_E b s (up 
to the energy reconstruction uncertainty of 25%). This 
conversion had been obtained as the average over 9 < 36° 
in Ref. [U using AIRES simulation code [|(J. That is, 
not only the average results were applied to individual 
showers, but effectively muon densities were simulated 
with CORSIKA while energies - with AIRES, though 
the two codes result in a systematically different rela- 
tions between energy and S^OO). Artificially high ener- 
gies resulted in higher, closer to observed, muon densities 
for simulated photonic showers. In our event-by event 
simulations with CORSIKA, the energies of gamma-ray 
primaries whose 5(600) were compatible to observed val- 
ues, were not higher by a factor 1.2, but in fact even lower 
than S b s for some of the events: besides the difference 
in simulation codes, this is partially due to non-uniform 
distribution of the highest-energy AGASA events on the 
celestial sphere ^3 US which makes the usage of aver- 
aged energies poorly motivated. 

The impact of two other sources of difference between 
our approach and that of Ref. [3(| is less important for 
the final result: (i) Ref. [3(| does not account for the 
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FIG. 4: Illustration of the difference between our study and 
Ref. |3(J. Plotted is the muon density at 1000 m versus the 
primary energy. Small grey boxes: simulated gamma-induced 
events with arrival direction of the event 1. Filled box, marked 
"simulated": simulated events compatible with the event 1 
by S^BOO). Open box, marked "RH": simulated showers with 
average Eo = 1.2_B b s from Ref. |30l| . The observed value of 
p,j(1000) for the event 1 is represented by a horizontal line, 
marked "observed". Eo is measured in eV, p,j(1000) in m -2 . 
See the text for more details. 

"lost" photons and (ii) the detector error is applied in 
our study to the simulated events while in Ref. [3(| - to 
the observed ones. 

The difference with Ref. |30j is illustrated in Fig. 0] 
where /^(lOOO) is plotted versus E for simulated 
gamma-induced showers with the arrival direction of the 
event #1. For simulated events compatible with the real 
event by S^OO), the average point is shown together with 
one sigma error bars. Horizontal error bars correspond to 
variations in Eq compatible with 5(600) . Vertical error 
bars include variations in simulated ^(1000) and 50% 
detector error. The point corresponding to simulated 
showers with Eq — 1.2E Q b s from Ref. [3(| has a larger 
p M (1000). Horizontal error bars correspond to the en- 
ergy reconstruction accuracy. Vertical error bars include 
variations in simulated /? M (1000) reported in Ref. [3(j and 
40% detector error applied to the observed value, added 



in quadrature. We see that the main source of the dis- 
agreement is in the values of Eq which push, for the case 
of Ref. [3(J, the simulated muon densities closer to the 
observed one. 



VI. CONCLUSIONS 

To summarize, we have studied the possibility that 
the highest-energy events observed by the AGASA and 
Yakutsk experiments were initiated by primary photons. 
Comparing the observed and simulated muon content of 
these showers, we reject this possibility for each of the 
ten events at E > 8 • 10 19 eV at least at the 95% C.L. An 
important ingredient in our study is the careful tracking 
of differences between the original and reconstructed en- 
ergies. This allows us to put an upper bound of 36% at 
95% C.L. on the fraction e 7 of primary photons with orig- 
inal energies Eo > 10 20 eV, assuming an isotropic photon 
flux and Eq 2 spectrum. This limit is the strongest one up 
to date. It strongly disfavors the Z-burst and constrains 
severely superheavy dark-matter models. The method 
that we have used is quite general and may be applied at 
other energies and to other observables. 
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